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We propose a method to decompose dynamical systems based on the idea that modules constrain 
the spread of perturbations. We find partitions of system variables that maximize ‘perturbation 
modularity’, defined as the autocovariance of coarse-grained perturbed trajectories. The measure 
effectively separates the fast intramodular from the slow intermodular dynamics of perturbation 
spreading (in this respect, it is a generalization of the ‘Markov stability’ method of network commu¬ 
nity detection). Our approach captures variation of modular organization across different system 
states, time scales, and in response to different kinds of perturbations: aspects of modularity which 
are all relevant to real-world dynamical systems. It offers a principled alternative to detecting 
communities in networks of statistical dependencies between system variables (e.g., ‘relevance net¬ 
works’ or ‘functional networks’). Using coupled logistic maps, we demonstrate that the method 
uncovers hierarchical modular organization planted in a system’s coupling matrix. Additionally, 
in homogeneously-coupled map lattices, it identifies the presence of self-organized modularity that 
depends on the initial state, dynamical parameters, and type of perturbations. Our approach offers 
a powerful tool for exploring the modular organization of complex dynamical systems. 


Many complex systems are modular, in that their com¬ 
ponents are organized in tightly-integrated subsystems 
that are weakly coupled to one another. Modularity has 
been argued to play many important roles, including in¬ 
creasing robustness CHa, evolvability ma, and func¬ 
tional differentiation [Si. Thus, there is great interest 
in measures of modularity and methods for decomposing 
complex systems into weakly coupled modules. 

This problem is here considered in the domain of mul¬ 
tivariate dynamics, a common formalism for modeling 
complex physical, biological, neural, and social systems. 
We propose a method of identifying dynamical modules 
motivated by the intuition that, in a modular system, 
the spread of perturbations is characterized by two time 
scales: fast spreading within modules and slow spreading 
between modules BIT]. In our treatment, the spread¬ 
ing process is coarse-grained relative to a partition (a de¬ 
composition of system variables into disjoint subsystems) 
by measuring the magnitude of the perturbation’s effect 
within each subsystem over time. If a partition reflects 
underlying modular structure, initially perturbed subsys¬ 
tems remain affected as dynamics unfold, while initially 
unperturbed subsystems remain largely unaffected. In 
this case, the partition’s coarse graining will capture the 
slow component of perturbation spreading dynamics, an 
effect quantified using a quality function called pertur¬ 
bation modularity. Our perturbation-based approach is 
related to many existing techniques for analyzing mul¬ 
tivariate dynamics, including Lyapunov-exponent based 
methods ISHIQ] and impulse response analysis m- 

As will be elaborated below, our methodology can 
identify the dependence of optimal decompositions on 
initial states, time scales, and kinds of perturbations ap¬ 
plied. These factors are all important aspects of mod- 
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ular organization in real-world dynamical systems. De¬ 
pendence on the initial condition reflects that dynami¬ 
cal systems can exhibit different modular organizations 
in different regions of their state-space; for example, 
distributed regions of the brain can couple into modu¬ 
lar assemblies via oscillatory synchronization, with the 
same region participating in different assemblies depend¬ 
ing on brain state [Tans]. The choice of time scale af¬ 
fects optimal decompositions by determining the sepa¬ 
ration between intramodular and intermodular pertur¬ 
bation spreading; in real-world complex systems, longer 
time scales have often been argued to correspond to 
larger-scale modules [DiiiHni- Finally, the dependence 
on the kinds of perturbations reflects that a dynamical 
system may be robust to some perturbations but highly- 
sensitive to others HU; for example, in biological double¬ 
knockout experiments, cellular responses to the simulta¬ 
neous deactivation of two genes can differ dramatically 
from responses to the individual deactivation of either 
gene [19]. 

Our approach starts from a pre-specified dynamical 
system and thus differs fundamentally from existing 
treatments of modularity based on network represen¬ 
tations of a system. Such methods are usually unable 
to capture the variation of modular organization across 
state-space or time scale, as well as other important dy¬ 
namical aspects of modularity. 

For instance, one standard approach applies graph- 
based community detection [20] to the structural network 
underlying a dynamical process (e.g., the social network 
over which an epidemic spreads). This treatment ignores 
the fact that the same structural network can support 
many different dynamical processes (for example, ‘com¬ 
plex contagion’ epidemics proceed differently from ‘sim¬ 
ple contagion’ epidemics [21]). In contrast, our method¬ 
ology is by definition sensitive to dynamical differences. 

Another class of methods detects community struc- 
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ture in network representations of dynamics, defined ei¬ 
ther in terms of causal interactions or statistical depen¬ 
dencies between variables (e.g., relevance networks in 
systems biology [22] and functional networks in neuro¬ 
science [23|). Unfortunately, constructing such networks 
involves a conversion of the dynamical system (defined 
in terms of transitions between multidimensional states) 
into a graph (defined in terms of nodes and edges). This 
conversion can affect modular decompositions in opaque 
ways as well as invalidate presumed graph-theoretic null 
models [24l |25] ; statistical dependency networks in par¬ 
ticular require a number of nontrivial decisions regarding 
the choice of dependency measure (correlation, transfer 
entropy, phase-locking measures, etc.), treatment of pos¬ 
itive versus negative interactions, and thresholding [23] . 
Furthermore, coupling between variables does not neces¬ 
sarily give rise to large values of correlation or other de¬ 
pendency measures [26| (also as shown in our first exam¬ 
ple below). Finally, community detection on dependency 
networks optimizes quality functions that are difficult to 
interpret in terms of the original system dynamics. Per¬ 
turbation modularity does not require the construction 
of a network representation of a dynamical system and is 
interpretable in terms of the separation of slow and fast 
time scales of perturbation spreading. 

Because our approach is based on intrinsic system dy¬ 
namics, it also differs from methods that identify modules 
by imposing a dynamical process onto a given network, 
such as diffusion of random walkers EZHU] or coupled 
phase oscillators [I5l|29]. However, as we discuss below, 
in certain cases our approach has connections to such 
methods. In particular, it can be seen as a generalization 
of the random-walk-based approach of Markov stability 
[28l l3Q] to a broad class of dynamics. 

To formally define perturbation modularity, consider a 
dynamical system with an A^-dimensional state-space X 
and evolution operator /^ : X ^ X at time scale t (both 
state and time can be continuous or discrete). Given a 
set £ of possible initial perturbations, s G f is applied 
to an initial condition cc G X to produce a perturbed 
initial condition x e G X. After time t, the size of 
the perturbation in the whole system is measured as the 
norm of the difference between the perturbed and unper¬ 
turbed trajectories: \\f^{x ^ e) — f^{x)\\. The relative 
size of the perturbation within a subsystem S (a subset 
of system variables) is: 


ms{x,£) 


\\(f(x + e)-fix))s\\ 

\\f{x + £) - p{x)\\ 


( 1 ) 


where the subscript S on the right hand side indicates a 
projection onto the dimensions indexed by S. For sim¬ 
plicity, we consider only cases where the system’s per¬ 
turbed and unperturbed trajectories have not merged 
{\\f^{x e) — f^{x)\\ > 0) and Eq. Q is well-defined. 

Assume a partition of the system tt = {Si,... ,Sk} 
into K disjoint subsystems. The coarse-grained perturba¬ 
tion vector yl^{x,s) = (cc, s),..., (cc, s)) cap¬ 
tures the relative size of the perturbation in each sub¬ 


system. Due to the normalization in Eq. vl is in- 
variant to the dynamical expansion of the whole-system 
phase-space, instead reflecting only the relative effects of 
perturbations on different subsystems. 

We now define perturbation modularity Q^(7r, x) as the 
vector autocovariance of the coarse-grained perturbation 
vector: 

Q*{tt,x)= ( 2 ) 

E[y°(®,e) •y*(a;,£)] - E[y°{x,£)] ■E[yl{x,e)] 

where the expected values are taken over P(s), a prob¬ 
ability distribution over perturbations (i.e. the elements 
oi £). The first term of Eq. 0 measures the degree to 
which perturbations persist within a partition’s subsys¬ 
tems (i.e. initially perturbed subsystems remain affected 
after time t, while initially unperturbed subsystems re¬ 
main relatively unaffected). The second term of Eq. 0 
provides a baseline expectation of perturbation effects 
that accounts for differences in subsystem sizes. 

As stated, the spread of perturbations in a mod¬ 
ular system will be constrained by module bound¬ 
aries. The optimal modular decomposition is the 
partition that maximizes perturbation modularity: 
M = arg max^^ (tt, cc). 

Perturbation modularity (Eq. 0 ), as well as opti¬ 
mal modular decompositions, are state-dependent in that 
they are defined relative to an initial condition x. Differ¬ 
ent criteria may be used to determine the choice of this 
initial condition, such as dynamical importance (e.g., an 
equilibrium state), particular research interest, or ran¬ 
dom selection. Alternatively, the modularity of entire 
state-space regions, rather than individual states, can be 
measured as the expectation of perturbation modularity 
over a distribution of initial conditions (e.g., by averaging 
across the entire system state space). Similarly, stochas¬ 
tic dynamical systems can be accommodated by taking 
expectations over future state distributions. Eor simplic¬ 
ity, however, these extensions are not considered in the 
present work. 

In addition to initial condition, perturbation modu¬ 
larity and optimal decompositions also depend on the 
time scale t, which, as mentioned, can act as a resolu¬ 
tion parameter. When there is not a time scale of a pri¬ 
ori interest, optimal decompositions can be identified at 
multiple resolutions by sweeping across a range of time 
scales. Einally, the measure also depends on the kinds 
of perturbations applied, as specified by £ and P{s). In 
practice, perturbations can be selected according to do¬ 
main knowledge (e.g., typically encountered environmen¬ 
tal perturbations) or using ‘neutral’ options (e.g., small 
increments to single variables). In many cases, initial 
perturbations should be localized to a small number of 
variables (i.e. the elements of £ are sparse) because the 
spread of perturbations is more pronounced when only a 
few subsystems are initially perturbed. As we will show, 
perturbations that simultaneously affect many variables 
probe the system at larger scales and uncover larger mod- 
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ules, providing another way to explore decompositions at 
multiple resolutions. 

Like other temporally-localized methods m, pertur¬ 
bation modularity also depends on the norm used to mea¬ 
sure perturbation magnitude. Below, the £i norm is used 
because it performs well and permits connections to com¬ 
munity detection methods in graphs (see Supplemental 
Material [32] for definition of £p norms). 

Perturbation modularity is related to the Markov sta¬ 
bility method of community detection in graphs, which 
identifies communities as subgraphs that trap random 
walkers [SHiisn]- Similarly to perturbation modularity, 
Markov stability separates diffusion dynamics into fast 
intracommunity and slow intercommunity components. 
As shown in the Supplemental Material [32], perturba¬ 
tion modularity is equivalent to Markov stability when 
the system of interest exhibits diffusion dynamics, per¬ 
turbations are homogenous increases to single variables, 
and the £i norm is used to measure perturbation mag¬ 
nitude. More broadly, our approach can be seen as a 
generalization of Markov Stability to other dynamics. 

In addition, £i perturbation modularity on a dynam¬ 
ical system is equivalent to directed weighted Newman’s 
modularity [33] [34] on a specially constructed graph (see 
Supplemental Material [32]). In this graph, nodes corre¬ 
spond to system variables and the edge from node i to 
node j has weight: 


Wi 


= E 


mliy{x,e) m\jy{x,e) 


where the expectation is over P{s) and the subscripts {i} 
and {j} indicate single-variable subsystems. This map¬ 
ping permits perturbation modularity to be maximized 
using highly efficient existing community detection algo¬ 
rithms (such as the Louvain method [35] [36] used for the 
examples below; code is available online m- 

Several criteria can be used to measure the quality of 
identified decompositions. High-quality decompositions 
have large perturbation modularity values (e.g., near 1 
for £i or £2 perturbation modularity, see Supplemental 
Material [32] for derivation of bounds on perturbation 
modularity). Additionally, high-quality decompositions 
are robust to small changes in system and optimization 
parameters. This can be quantified by measures of parti¬ 
tion similarity like normalized mutual information (NMI) 
[38] . an information-theoretic measure that ranges from 
0 (maximally dissimilar partitions) to 1 (identical parti¬ 
tions). In several of the examples below, we plot NMI 
similarity between optimal decompositions identified at 
close values of t; high NMI values indicate modular or¬ 
ganization robust to small changes in time scale. Similar 
techniques are used in the Markov stability literature to 
identify time scales with robust decompositions [39] . 

We demonstrate our method on several examples of 
coupled logistic maps, nonlinear discrete-time dynami¬ 
cal systems that have been used to explore spatially- 
extended chaos and pattern formation [40]. Assume a 
system of N variables, with Xi{t) indicating the state of 
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Figure 1. (Color online) System of 80 coupled logistic maps 
(0=2, 7=0.04) (a) The coupling matrix exhibits hierarchical 
modularity at three levels, (b) The system is chaotic and 
no strong correlations between variables are observed over 
10,000 time steps, (c) Perturbation modularity of optimal 
decompositions (PM, solid black) at different time scales t 
and NMI between optimal decompositions at successive times 
(dashed blue). Three time scale regions are robust (NMI=1, 
grey), corresponding to the three hierarchical levels of the 
coupling matrix (insets). 


variable i at time t, and the transition function: 

Xi{t + 1) = {1 - (3) 

where g{x) = 1 — ax‘^ is the logistic map, parameter 
a e [ 1 , 2 ] controls the chaoticity, parameter 7 G [ 0 , 1 ] con¬ 
trols the coupling strength, ‘coupling matrix’ elements 
kji determine the influence of variable j on variable i, and 
di = kji normalizes the coupling strengths. When 

variables are homogeneously coupled to nearest neigh¬ 
bors on a 1 -dimensional ring lattice, these systems are 
called coupled map lattices (CML) [40]. Coupled logistic 
maps display a rich variety of spatiotemporal patterns 
in different parameter regimes due to the interplay be¬ 
tween intervariable coupling (which ‘homogenizes’ vari¬ 
able states) and chaos (which injects variation into vari¬ 
able states). 

We consider several examples of coupled logistic maps. 
Unless otherwise noted, perturbations consist of a uni¬ 
form distribution over small increases to single vari¬ 
ables: £ = { 0.0001 ’ Bi : i = l.W}, where is the i^^ 
Wdimensional standard basis vector. The £i norm is 
used to measure perturbation size. 

In Example we uncover modular organization that 
is present in a system’s coupling matrix, though not ap¬ 
parent in the correlation statistics. Consider an ^"=80 
variable system with chaotic dynamics ((a=2,7=0.04) 
and a hierarchical modular coupling matrix (Fig. El)- 
The system is composed of 8 tightly-coupled low-level 
modules {kji=l) with 10 variables each, pairs of which 
are nested within 4 mid-level modules {kji=0.01), pairs 
of which are in turn nested within 2 weakly coupled high- 
level modules (/cj 4 = 0 . 0001 ). A random state is used as 
the initial condition. 

Because the system is strongly chaotic for these val- 
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Figure 2. (Color online) Two 100-variable CMLs are com¬ 
pared: one ‘modular’ (top row; a=1.7,7=0.1) and one ‘diffu¬ 
sive’ (bottom row; «= 1.9,7=0.6). (a,d) Spacetime plots of 
the effect of a single-variable perturbation. A pixel is colored 
black if the absolute difference between perturbed and unper¬ 
turbed trajectory at a variable (vertical axis) exceeds 1% of 
the size of the system-wide perturbation at a given time (hor¬ 
izontal axis). (b,e) Spacetime plots of the optimal decompo¬ 
sitions at different time scales; color indicates each variable’s 
subsystem. (c,f) Perturbation modularity of optimal decom¬ 
positions (PM, solid black) at different time scales t and NMI 
between optimal decompositions at successive times (dashed 
blue). Stable decompositions are observed in the modular 
CML (top row). 


ues of a and 7 , there is no obvious ‘order parameter’ 
for identifying modular organization from system trajec¬ 
tories US]; for instance, variable states are largely un¬ 
correlated over 10,000 time steps (Fig. Ef)- However, be¬ 
cause perturbations first spread within low-level modules, 
then mid-level modules, and finally high-level modules, 
our method easily uncovers the hierarchical modular or¬ 
ganization. Fig. shows the perturbation modularity 
(black) and NMI (dashed blue) for optimal decomposi¬ 
tions at different time scales. There are three robust time 
scale regions, corresponding to each of the three hierar¬ 
chical levels of the coupling matrix (insets in Fig.[^). Be¬ 
yond time scale ^50, perturbations have spread between 
the high-level modules; at this point, optimal decompo¬ 
sitions reflect random fluctuations in initial conditions, 
and perturbation modularity and NMI values are near 0. 

In Example we investigate a more interesting case 
in which modularity emerges in a homogeneously-coupled 
CML. In some parameter regimes, spatial variation in ini¬ 
tial conditions breaks the lattice coupling symmetry and 
leads to the emergence of modular domains (contiguous 
lattice regions) that constrain the spread of perturba¬ 
tions [41]. Such a ‘modular’ regime is investigated using a 
CML with N=100 variables and weak coupling-strength 
((a=1.7, 7 = 0 . 1 ). The initial condition is set by iterating 
a random state for 10,000 time steps. Fig. shows the 
spacetime plot of the effect of a single-variable perturba¬ 
tion to this initial condition: the perturbation spreads to 
several nearby variables until t^50 but then remains con¬ 
fined within its domain. When computed over a uniform 
distribution of single-variable perturbations, our method 
uncovers robust modular organization for a large range 
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Figure 3. (Color online) State-dependence in the modular 
CML. For a 12,000 step trajectory starting from a random 
state, optimal decompositions (time scale t=300) are iden¬ 
tified using states along this trajectory as initial conditions, 
(a) Optimal perturbation modularity (PM, solid black) grows 
with increasing trajectory steps (horizontal axis), indicating 
the emergence of robust modular structures. Trajectory step 
10, 000 (dotted green line) is the initial condition in examples 
[^ and (b) Optimal decompositions identified at different 
trajectory steps; color indicates the subsystem of each vari¬ 
able (vertical axis). 


of time scales (Fig. If)’ with optimal decompositions ex¬ 
hibiting high values of perturbation modularity and NMI 

(Fig#). 

The above system can be compared to a CML in a ‘dif¬ 
fusive’ regime ((a=1.9, 7 = 0 . 6 ). For these parameters, the 
effects of perturbations spread freely across the lattice, as 
shown in the spacetime plot of Fig. (initial condition 
is the same random state as in the modular CML iterated 
for 10,000 time steps). This system does not exhibit ro¬ 
bust modular organization: optimal decompositions are 
not stable (Fig.[^) and optimal perturbation modularity 
and NMI values are low (Fig.|^). Once the effects of per¬ 
turbations spread completely around the ring lattice at 
t^lOO, both optimal perturbation modularity and NMI 
values are near 0. 

In Example we demonstrate state-dependent mod¬ 
ularity by tracking the gradual emergence of modular 
organization over the course of a CML trajectory. The 
modular CML of Example ((a=1.7,7=0.1) was started 

from a random state and iterated over a 12, 000 step tra¬ 
jectory. The state encountered after 10, 000 time steps 
was previously used as the initial condition in Example 
Here we find optimal decompositions (time scale t=300) 
when different states along the aforementioned trajectory 
are used as initial conditions. Over the course of the tra¬ 
jectory, optimal perturbation modularity grows through 
a series of plateaus (Eig.[^), indicating the appearance of 
modular structures. EigT]^ shows the optimal decompo¬ 
sitions identified at different trajectory steps. Variables 
^1—40 quickly form modular structures (by trajectory 
step ^2000), while variables ^40—100 need almost 7, 000 
steps to do so. This provides an example of self-organized 
modularity^ or modular organization that emerges during 
a system’s dynamical evolution. 

Previously, we showed that perturbation modularity 
captures the presence of stable modular structures in dif¬ 
ferent CML parameter regimes (Example]^, and that it 
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Figure 4. (Color online) Parameter phase map of CML. Per¬ 
turbation modularity (color) for optimal decompositions of 
100-variable CMLs with different values of chaoticity (a; hor¬ 
izontal axes) and coupling strength ( 7 ; vertical axes). Pertur¬ 
bation modularity is computed at three time scales for two dif¬ 
ferent classes of initial conditions: random states [(a) t= 100 , 
(b) t=200, and (c) t=300] as well as random states iterated 
for 10,000 time steps [(d) t=100, (e) t=200, and (f) t=300]. 


can uncover modular organization in a state-dependent 
manner (Example [^. In Example we use pertur¬ 
bation modularity to characterize regions in the CML 
parameter space with respect to the onset of modularity. 

Specifically, we construct 100 -variable CMLs with dif¬ 
ferent values of chaoticity (a) and coupling ( 7 ) parame¬ 
ters. For these different parameter values, Fig. shows 
optimal perturbation modularity computed at three time 
scales (t=100, t=200, and t=300) and two different 
classes of initial conditions: random states (Fig. l^-c) 
and random states iterated for 10,000 time steps ^ig. 

-f). In all cases, optimal perturbation modularity val¬ 
ues were averaged across 10 random samples of initial 
conditions. 

Several regimes of spatial organization can be iden¬ 
tified in the parameter phase maps. For a < 1.44, 
spatial domains, which form even when the system is 
started from random initial conditions, constrain the 
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Figure 5. (Color online) Perturbation dependence in the mod¬ 
ular CML. (a) Spacetime plot of the effect of a perturbation 
to a 20-variable window (red arrows), as in Fig. and |^. 
(b) Optimal decompositions (time scale t=300) for different 
perturbation sizes (horizontal axis). 


spread of perturbations over long time scales; we call 
this the modular regime. For other parameter values 
(e.g., 1-6 < Of < 1.95, 7 ^ 0.1, the yellow ‘tongue’ in 
Fig. ^ -f), modular domains only appear when random 
states are iterated for many steps before being used as 
initial conditions. This regime, which includes the case 
studied in Example we call self-organized modular. Fi¬ 
nally, for parameter values corresponding to the blue re¬ 
gions in Fig. which we call the diffuse regime, mod¬ 
ular domains are not present and perturbations spread 
freely. Here, different parameter values give rise to dif¬ 
ferent diffusion speeds [4Tj: for example, (a= 1 . 9 , 7 = 0.7 
exhibits no modular organization at time scale t= 100 ; 
on the other hand, 0 ^= 1 .9, 7 = 0.2 maintains some mod¬ 
ularity at t= 100 , but this organization disintegrates by 
t=200. 

Finally, in Example we explore modularity’s de¬ 
pendence on the kinds of perturbations applied. We 
again consider the modular CML ( 0 ^= 1 .7, 7 = 0 . 1 ) and 
initial condition of Example Instead of perturbing 
single variables, we now simultaneously perturb multi¬ 
ple variables in lattice-contiguous ‘windows’ of different 
sizes (variables simultaneously incremented by 0 . 0001 ; all 
N windows are perturbed with uniform probability); for 
illustration. Fig. lit shows the effect of a perturbation to 
a window of 20 variables. Fig. shows that optimal 
decompositions (time scale t=300) depend on perturba¬ 
tion size. As more variables are perturbed, smaller sub¬ 
systems merge into larger subsystems in a hierarchical 
manner. 

Future work can pursue several extensions to our ap¬ 
proach. First, estimating perturbation modularity from 
real-world datasets is of great practical interest; this can 
be investigated by applying the method to fitted dy¬ 
namical models (e.g., vector autoregressive or dynam¬ 
ical causal modeling [42]) or using nonparametric ap¬ 
proaches. Second, it is possible to explore other mea¬ 
sures of decomposition quality beyond robustness to time 
scale, including robustness to changes in initial condi¬ 
tions and kinds of perturbations; alternatively, decompo¬ 
sition quality may be evaluated by testing the statistical 
significance of optimal perturbation modularity against 
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null-model ensembles of nonmodular dynamical systems 

[24] . Third, it is of interest to identify possible limita¬ 
tions of our method, such as whether it is susceptible 
to the kinds of resolution limits [43] and detectability 
thresholds [44] encountered by graph-based community 
detection methods. Finally, future research can investi¬ 
gate other measures of perturbation magnitude (e.g., dif¬ 
ferent norms or divergence functions), kinds of decompo¬ 
sitions (e.g., overlapping subsystems), and cost functions 
(beyond vector autocovariance). For example, cost func¬ 
tions that capture the invertibility or sparsity of coarse¬ 
grained dynamics could be used to decompose a system 
into a mesoscopic ‘control diagram’, in which each sub¬ 
system controls a small number of others. 

To summarize, we identify modular decompositions of 
multivariate dynamical systems based on the idea that 
modules constrain the spread of perturbations. We pro¬ 
pose a quality function, called perturbation modularity. 


which can be used to identify optimal coarse grainings 
that capture the slow component of perturbation spread¬ 
ing dynamics. The method generalizes graph-based com¬ 
munity detection to a broad class of nonlinear dynamical 
systems and provides a principled alternative to detect¬ 
ing communities in network representations of dynamics. 
The method captures variation in modular organization 
across different time scales, initial conditions, and kinds 
of perturbations and offers a powerful tool for exploring 
modularity in complex systems. 
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Supplemental Material 

Appendix A: Bounds on Perturbation Modularity 


The bounds on perturbation modularity depend on the norm used to measure perturbation magnitude (see Eq. (1) 
in the main text). Here we consider the family of ip norms, where the ip norm of vector a = (ai,..., uat) is defined 

for p > 0 as ||a||p = |a/) ' . 

As we show below, perturbation modularity is bounded by —1 and 1 for ii and ^2 norms. More generally, for ip 
norms withp > 2, perturbation modularity is bounded by and where N is the number of variables 

in the original system. 

Assume some ip norm is used to measure perturbation magnitude. Using the definition of perturbation modularity 
for initial condition x and partition tt: 

Q\Tr,x) =E[y°{x,£) ■yl{x,£)] -E[y°(a;,e)] ■E[y*(a;,e)] 

<E[y°{x,£) ■yi{x,£)] 

<E[||y°(a:,e)||2||y*(a;,e)||2] 

where the second line follows from the non-negativity of the coarse-grained perturbation vectors y^{x, s) and yl(x, s) 
and the third line follows from the Cauchy-Schwarz inequality. 

When some ip norm is used to measure perturbation magnitude, the coarse-grained perturbation vectors themselves 
are unit vectors in ip. To show this, let Vi = l(f^(x -he) — |. Then, 




VS'Gtt 


V 1/p 

/ ■ 

= 

E 

/ 

\^se7T 


-1 px 1/p 


(E, 


es 


"fV 


.{Y.L,4y. 


= E& 


1/p 


es ' 


KSen 




= 1 


When p = 2 is used to measure perturbation magnitudes, the upper bound can be rewritten as: 

g*(7r,a:) < E[||y°(a;, e )||2 ||y* (a;,e)|| 2 ] = E[1 ■ 1] = 1 


When p = 1 is used to measure perturbation magnitudes, we first note that ||a ||2 < ||ot||i for any a and that 
p^(cc,s)||^ = ||p^(cc, s) 11^ = 1. The upper bound can be rewritten as: 

g*(7r,a:) < E[||y°(a;,e )||2 ||y^(a;,e)|| 2 ] < E[1 ■ 1] = 1 


More generally, whenp > 2, Holder’s inequality provides the bound ||a ||2 p) l|o||p) where n is the number 

of dimensions of a. Thus, ||p^(cc, s )||2 < kl*'"' ^i\yi(.x,£)\\j^= p)<Ar(2 p), since N is the maximum number 

of subsets in a partition. This gives: 


g*( 7 r,a;) < E[||y°(®,£)||2 ||y* (®,£)||2] < E 


iv(5“p)77(5-^)l = N^er) 


To show the lower bound, we again use the definition of perturbation modularity and similar reasoning: 

Q\'k,x) =E[y°{x,£) ■yi{x,£)] -E[y°(a;,£)] •E[yka;,£)] 

> -E[y°(®,£)] ■E[yl{x,£)] 

> - || E [ y °(«, e )]||2 || E [ yka ;,£)]||2 

> -E[||yO(a::,£)||2]E[||y°(a;,£)||2] 

where the last line follows from Jensen’s inequality and the convexity of norms. Similar arguments as above show 
that (tt, cc) > — 1 for p = 1 and p = 2 and (tt, x) > for p > 2. 

In practice, we are also interested in the maximal perturbation modularity across all partitions. In fact, maximal 
perturbation modularity is always non-negative. This is because there is at least one partition with perturbation 
modularity equal to 0: the partition that contains the entire system as one subsystem. In this partition, which we 
call TTo = {{1,..., N}}, perturbations are entirely contained in the single subsystem and y^^{x^ s) = y'^^{x, s) = (1) 
for all X and s. From the definition of perturbation modularity, it can be seen that (tto, x) = 0 for all t. 









Appendix B: Perturbation modularity and community detection 


An explicit connection can be made between our approach and graph-based community detection. We first rewrite 
and expand the definition of perturbation modularity from the main text: 

Q*(7r,a;) = E[t/°(a:,e) ■ y* (a:,e)] - E[y°{x,£)] •E[y^(a;,e)] 

= ^ {E[m%{x,£)m%{x,e)] - E[m^(®,£)] E [m%{x,£)]) 

SEtt 

where the expectations are over P{s). We assume that the £i norm is used to measure perturbation magnitude, which 
provides the following additive property: mg{x,s) = where the subscript {i} indicates a subsystem 

only containing variable i. We rewrite the above equation as: 


Q\7T,x) = Y^ E 


Sen 


m\^y{x,e) 


ies 


jes 


-E 




ies 


E 




jes 




Sen \i,jes ies jes 

where the last line comes from exchanging the order of summation and expectation. 


(Bl) 


Bl. Perturbation modularity is Markov Stability for diffusion dynamics 


Markov stability is a method of community detection that uses the dynamics random walks over graphs [SHlEol. 
Here, communities are defined as subgraphs that tend to trap random walkers. This method is of particular interest 
because it generalizes many other community detection methods, including optimization of Newman’s modularity, 
cut size, and spectral methods [30]. Given a random walk over an A-node graph, the Markov stability of a partition 
TT at time scale t is defined as: 

H'M = E Pr(walker in S at times 0 and t) — Pr(walker in S' at time 0) Pr(walker in S at time t) (B2) 

Sen 

In this framework, the optimal partition of a graph is the one that maximizes Markov stability. 

Given homogenous perturbation to single variables, there is an equivalence between Markov stability and £i per¬ 
turbation modularity of diffusion dynamics. Specifically, the diffusion of random walkers on a graph can be stated in 
terms of a linear dynamical system (cc) = T^cc, where x (t) is an A-dimensional vector with Xi (t) being the density 
of random walkers at node i at time t, and is the time scale t transition matrix (here superscript t indicates matrix 
power; for simplicity, we consider the discrete-time case). Assume that perturbation to variable i is indicated by 
s* = e^, where is the i^^ A-dimensional standard basis vector (initial perturbations can be scaled by any constant 


{j}| _ _ X. 

^^Ili “ lledli 

is a transition matrix: it has positive entries and preserves the £i norm upon matrix multiplication. Therefore: 


without changing results). Then, m|^.|(cc, s*) = 
rtrix: it has positive < 

(/*(*+ 60 -/*(*)){,.} 




\{T\x + e^)-T^x){^}\ \{Te%y\ 


\\P{X + £^) - fix) 


\\Tt{x + £i)-Ttx\\^ 


\\Tte% 




_ ij^t 




The terms in Eq. (Bl) can now be mapped to the terms in the Markov stability Eq. 

N 

P{e'^)Tk = Pr(walker in S at times 0 and t) 

i,jeS iJ^S k=l iJ^S 

N 


E^t 

ies 




=EE P{£'')Si,k = E P{s^) = Pr(walker in S at time 0) 

ies k=i ies 


N 


Ee rn\j}{x,£) = P{£^)Tl.j = Pr(walker in S at time t) 

jes jeSk=i 
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Markov stability is usually defined for an equilibrium random walk, such that p (walker in S at time 0) = 
p (walker in S at time t). In our framework, this is accomplished by setting P (s*) equal to the equilibrium probability 
of finding a random walker at node i. 


B 2 . Mapping to directed weighted Newman’s modularity 


In this section, we show that for any dynamical system, ii perturbation modularity can be mapped to directed 
weighted Newman^s modularity on a specially-constructed graph. One result of this mapping is that efficient community 
detection algorithms can be used to find partitions that maximize perturbation modularity. 

To review, the weighted directed Newman’s modularity of partition tt is defined as [33]: 


E E 

C GTriJ eC 




M 


(B3) 


where Wij indicates edge weight from node i to node j wij is the out-degree, Wij is the in-degree. 


and M = Wij = is the total strength (summations in these equations are over all nodes). 

Now, for an A^-dimensional dynamical system, we construct a graph with one node for each variable and edge 
weight from node i to j: 


Wij = E 




When the £i norm is used, = 1- Thus, 
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Rewriting Eq. (Bl) makes the mapping to Eq. (B3) explicit: 


Q*(7r,®) = E E ^ mliy{x,£)m\jjx,£) -^^E m\ijx,£)^ ^^E m\jjx,£) 
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